Hydrodynamic Tensor-DFT with correct susceptibility 
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In a previous work we developed a family of orbital- free tensor equations for DFT [J. Chem. 
Phys. 124, 024105 (2006)]. The theory is a combination of the coupled hydrodynamic moment 
equations hierarchy with a cumulant truncation of the one-body electron density matrix. A basic 
ingredient in the theory is how to truncate the series of equation of motion for the moments. In 
the original work we assumed that the cumulants vanish above a certain order (N). Here we show 
how to modify this assumption to obtain the correct susceptibilities. This is done for N=3, a 
level above the previous study. At the desired truncation level a few relevant terms are added, 
which, with the right combination of coefficients, lead to excellent agreement with the Kohn-Sham 
Lindhard susceptibilities for an uninteracting system. The approach is also powerful away from 
linear response, as demonstrated in a non-perturbative study of a jellium with a repulsive core, 
where excellent matching with Kohn-Sham simulations is obtained while the Thomas Fermi and 
von-Weiszacker methods show significant deviations. In addition, time-dependent linear response 
studies at the new N=3 level demonstrate our previous assertion that as the order of the theory is 
increased, new additional transverse sound modes appear mimicking the RPA transverse dispersion 
region. 
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I. INTRODUCTION 

The development of new methods for quantum dy- 
namics based upon hydrodynamic representations is very 
promising. In hydrodynamics the kinetics of the system 
is defined by a lesser number of variables than the number 
of variables required to define the complete one-particle 
density matrix (which contains all the information on off- 
diagonal quantum coherence as in, e.g., the Kohn-Sham 
approach) . For stationary studies the hydrodynamics ap- 
proach is related to orbital-free density-functional theory 
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0J, Ha, [HI, [H, [H, 1H HI- :t is the reduced number of 
variables depicting the system that makes hydrodynami- 
cal theories applicable for numerical studies of relatively 
large systems. 

The simplest hydrodynamical approach is similar to 
the de Broglie-Bohm formulation of one-particle quan- 
tum mechanics [H M, M, HI M EH El El- In this 
approximation the complete complex-valued one-particle 
density matrix is substituted by two real valued fields 
p and <j), which are combined in an order parameter 
ip = ^fpexp(i(p). The equations of motion are obtained 
by minimizing a Ginzburg-Landau-like functional on tjj. 
In addition the density matrix is assumed to possess long- 
range off-diagonal one-particle correlations. 

A more rigorous and asymptotically exact approach is 
an infinite hierarchy of coupled hydrodynamic moment 
(CHM) equations [3, EE EE El El- The moments 
come from a Taylor expansion of the one-particle den- 
sity matrix with respect to the off-diagonal variable. To 
get a tractable system of equations the infinite hierar- 
chy must be truncated. The most physically meaningful 
truncation is a cumulant expansion for the density matrix 
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|34j . Specifically, one decides on an order to terminate 
the method at; a low order will be less numerically de- 
manding but less accurate than a higher one. Then, at 
that order, labeled N, the (N+l)-th order moment is ex- 
panded in terms of the previous set of moments, through 
the use of the cumulant expansion. 

The CHM theory and the accompanying cumulant 
truncation have been applied so far to systems where 
particle statistics does not play an important role. In 
Ref. [H we have generalized the CHM theory and cumu- 
lant expansion to statistically degenerate fcrmions. The 
main point has been the modification of the unperturbed 
one-particle density matrix of a locally homogenous elec- 
tron gas by using the cumulants. Since the approach 
uses successive tensors, we labeled it Hydrodynamic ten- 
sor DFT (HTDFT). 

It turns out that the lowest level of truncation, N = 1, 
HTDFT corresponds to a de-Brogilie-Bohm quantum hy- 
drodynamics and in addition naturally incorporates the 
Thomas-Fermi 40] kinetic energy term into the energy 
functional. At the next level, N = 2, HTDFT starts re- 
producing the spectrum of a homogenous Fermi liquid, 
i.e., it gets transverse excitations, rather than just clas- 
sical plasmonic longitudinal excitations. The transverse 
sound mode mimics the elementary excitations' density 
of states. 

A crucial feature of HTDFT is the value of the cumu- 
lant used at the truncation. In Ref. \^ we assumed that 
the (N + l)-th order cumulant is zero. It turns out, how- 
ever, as we show here, that this assumption leads to a 
wrong susceptibility for a homogenous electron gas, i.e., 
to a wrong linear response to a perturbation, even for a 
non-interacting system of electrons. We show here how to 
remedy this problem. This is exemplified below for trun- 
cation at the N — 3 level, which is the first level where 
the method will yield different ground-state results from 
the Thomas-Fermi approach. Specifically, the 4'th or- 
der cumulant is written as a sum of terms involving the 
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gradients of the previous moments. The coefficients of 
these terms are obtained by fitting to the exact suscepti- 
bility of a non-interacting set of electrons (the Lindhard 
function). 

The balance of the paper is as follows. The general 
methodology is first developed in Section II. In Section III 
the derivation of a correct susceptibility is done. Section 
IV applies the methodology to a static non-perturbative 
numerical study of a jellium with a deep spherically sym- 
metric hole, where we show that the agreement with 
Kohn-Sham results is excellent while the Thomas Fermi 
and the von-Weiszacker methods have significant errors. 
Section V is a linear-response time-dependent study of 
the approach for N=3 as a function of frequency and 
wavevector. This latter part is a direct continuation of 
our work in Ref.[39j for N = 2, and proves that there 
is an additional sound mode with respect to the N = 2 
case, just as suggested in Ref.[3^|. Conclusions follow in 
Section VI. 



external potential, and E xc is the exchange-correlation 
energy and po is the nuclei density. There are a variety 
of functions V xc = 5E XC /Sp in the literature (see, e.g., 
Ref. 0, H, 0]). For us, however, the specific form of V xc 
is not important. [In future works we will aim to derive 
a form of V xc which depends also on other moments in 
addition to p{R)] 

The particle kinetics in the system can be exactly de- 
scribed by the complete infinite set of hydrodynamic mo- 
ments (dynamic tensors) [13, HH, [H, H3, [H, HI] , which 
are the derivatives of the one-particle density matrix with 
respect to the off-diagonal distance, s, at s = 0: 



s=0 



(5) 



The particle and the current spatial densities are merely 
the first two tensors in the family: 



II. THE SYSTEM AND TENSOR-DFT 
FORMULATION 

A. Coupled Hydrodynamic Moment Hierarchy 



(6) 



By using Eq. © one derives an infinite set of equations 
which connects the moments at different orders: 



For completeness, we rederive the basic aspects of the 
theory (see Ref . 39] ) . We assume that the many electron 
system can be described by the one-particle density ma- 
trix, pW. The one-electron Hamiltonian governing this 
system, h, is, as usual, composed of kinetic terms, and a 
local potential terms. 

The one-particle density matrix is then expressed in 
terms of average and difference coordinates as: 

pW(R,s) = (ft(R-s/2)$(R + a/2)). (I) 

The time evolution of the one-particle density matrix is 
governed by the Heisenberg equation, ip = [h,p], which 
in those coordinates takes the form: 

at 

+ (V(R + a/2) - V(R - a/2)) P [1 \2) 

Here P a and p a stand for the derivatives over the coor- 
dinates R a and s a 



P a = —id/dR a , p a — —id/ds a , 



(3) 



and V(R) is the effective potential, which also takes into 
account the two-body interactions: 



v(R) = 



p(R')-po(R') d3Rf _ SE XC 



R R' 



Sp(R) 



■^(fl),(4) 



where p(R) = p^'(R,0) is the spatial electron density, 
Po(R) is the positive nuclear charge density V ex t is any 
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etc. 







(7a) 
(7b) 



(7d) 



This generic set of equations is correct for both fcrmions 
and bosons. For this set to be useful one should terminate 
it at some level. As usual, this termination is actually a 
method for factorizing a moment $( iv+1 ) at some N into 
moments fc < N. In addition, this truncation reflects 
the Fermi statistics of the particles. The order N at 
which one terminates controls the precision with which 
we treat the system. 



Fermi-factorization of higher order dynamic 
tensors 



In Ref.[39j we proposed a factorization procedure for 
the lowest order dynamic tensors (N — 2, 3). Here we de- 
scribe in detail how the factorization of the higher order 
dynamic tensors works in the Fermi case. The method 
proposed is based on the following general parametriza- 
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tion of the one-particle density matrix: 
p w (R,s) = peM<t>(R,s)}f (p,s), (8a) 
0(B,s) = ^^4>£l.. i jR)(is i J(is i2 )...(isM 



fo(p,s) 



a>l 

sin(fci?s) — (hps) cos(kps) 



(k F s) 3 



k F = (S^p) 1 / 3 . 



(8c) 
(8d) 



Here, /o is the normalized one-particle density matrix of 
a free fermion liquid with density p(R) and hp is the 
local (density-dependent) Fermi wave-vector. All the cu- 
mulants, 4>^ a \ are symmetric in all the indices because 
they are convolved with the symmetric tensors Si t . . . Si n 
and all the ( - a - ) are real as the one-particle density matrix 
is hermitian. The same is true for the tensors 



The physical meaning of this paramerization is as fol- 
lowing. If = 0, then we end up with the Thomas-Fermi 
approximation of a locally homogeneous Fermi liquid. 
The function perturbs this steady liquid picture, and 
the tensors Q 's and/or <j)( Q )'s determine different dy- 
namic characteristics of the flowing electron liquid. The 
function /o assures the Fermi statistics of the particles at 
the one-particle level. 

For brevity we introduce below the tensors 



(9) 



instead of $ (q) . The tensors and (a) are interre- 
lated. The relations between ^"w's and 0( a )'s for the 
lowest order tensors are given below for the first four re- 
lations: 



= (1) , (10a) 

jr( 2 ) = 0( 2 ) + 0(i)0(i) +e (2 \ (10b) 

jrO) _ 00) + 30(1)0(2) + 3 0(i) e (2) + 0(i)0(i)0(i) ; (10c) 

;r( 4 ) = 0W + 40(!)0(3) + 30 (2) (2) + 60( 1 )0( 1 )0(2) 



+60( 2 )e( 2 ) + 60( 1 )0( 1 )e< 2 ) + 0( 1 )0«0( 1 )0( 1 ) + e (4) , (lOd) 

etc. Here all terms are absolutely symmetric tensors of their indices so that there is no need to write down the indices 
explicitly; a bar denotes here a complete symmctrization, e.g., for a product of a^...^ and bi 1 ___i M : 

abi 1 ...i K+M — j jj^yj Q>p(i 1 )...p(i K )bp(i K+1 )... p (i K+M ), 

where summation is assumed over all (K + M)l permutations of the indices and p denotes a permutation. The 
symmetrized multiplication is associative and it can be considered a multiplication on a ring of symmetric tensors 
(note that in Ref. [391 ] a somewhat different symmetrization was used) . Here is an explicit example of the symmetrized 
multiplication: 

www - ~ (tMWf + tiWW + tfWffl 

In Eqs. (fTU| the tensors e' 2 ) and e^ come from differentiating the function /q, so that they have the physical meaning 
of averaging the particle momenta products over the unperturbed Fermi sea (ufs): 

4f = P~ % (PiPj)ufs = C 2%> (12) 



(4) _ -1 



jkl = P (PiPjPkPl) u f s = 3C4§Sijkl = C 4 (SijSkl + SikSjl + SilSjk) , (13) 

where the kinetic coefficients are defined as 

c 2 = \k%, C4 = ^4. (14) 

All the odd order e's vanish. 

The general recipe for how to express !F^ N ^ in terms of 0( Q )'s is as follows. is the sum of all the different 

symmetrized (in the sense discussed above) products of 0" and e^ a \a < N. An additional rule is that each term 
may include one (and only one) e-tensor. 
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The relations inverse to Eqs. (|T0|) are: 

JF« - 0« = 0, (15a) 

jr(2)_0(2) = jr(i)jr(i)+ e ( 2 ), (15b) 

J7(3)_0(3) = 3^(1)^(2) _ 2^(1)^(1)^(1) ; (15c ) 
jr(4) _ 0(4) = 4 jr(i)jr(3) _ 12^(1)^(1)^(2) + 6J r(i)jr(i) JT(i)jT(i) 

+3J C "( 2 )J C '( 2 ) -3e( 2 )e( 2 ) +e (4) . (15d) 



The inversion of the infinite set of relations (TlU|) is pos- 
sible since an expression for any in terms of 's 
contains only <f>^ a ', a < N. This means that if one 
knows the expressions for the first N cf)^ a > tensors {e.g., 
Eqs. (|15alll5cp ) for N — 3) then by substituting all the 
lower order 0( Q )'s, a < N, in the relation for J r ( N+1 ) 
(Eq. ljlOdjl ) with corresponding expressions in terms of 
s one gets the inverse relation for cf>( N+1 ) (Eq.(Tl5d|). 
The factorization for the tensor $( w+1 ) is then simply 
given by the (N + l) th equation in Eqs. [15]. The expres- 
sion for the $( w+1 )' tensor contains only kinetic tensors 
of order n < N , as well as the (N + l) th order cumulant. 
Once this cumulant is known the system of Eqs. (8) closes 
and one arrives at the N th order tensor DFT theory. 
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C. JV = 3 Hydrodynamic tensor-DFT 

At the N = 3 level the factorization of <i>( 4 ) is given by 
Eq. (|15d[) with (f)^ set to zero, or: 



$ (4) = 4p^ J$( 3 ) - 12p~ 2 JJ$( 2 ) 
+6p~ 3 JJJJ 



+3 i o- 1 $( 2 )$(2) - 3pe( 2 )e( 2 ) + peW + p0 (4 ^16) 

In order to complete the theory we need to obtain (j)( 4 \ 
For this, we study the static linear response of a homoge- 
neous electron gas. In the ground state all the odd order 
$( Q ' -tensors vanish (a ground state has no currents as its 
wave-function is real when there is no magnetic field and 
no degeneracy), so that the first three terms in Eq. (fT6|) 
would give only non-linear contributions and can be ne- 
glected. As a result, the required factorization for static 
studies simplifies as (here we restore the indices): 



(4) 
ijkl 



+3p(c4p)-4(p))S ij S k i+M l . (17) 



FIG. 1: The bare structure factors, Xo, on a scaled momentum 
scale, q/2kF, for the 1/9-von Weiszacker approach, bare and 
adjusted $ ( - 3) HTDFT theories, and the Lindhard function 
(free-electron gas static density-density correlator). The $' 3 ' 
HTDFT is fitted to have the three properties of the Lindhard 
function given in Eq. (|27|l . 



The structure factor is actually the static limit 
of the density-density correlation function x(q) — 
(p(-q, -w)/5(q,w))|^ i0 . The physical meaning of x(q) 
is the ratio between the amplitude of the infinitesimal 
harmonic change in electron density, p(q), and that of 
the external potential, v ext (q), which induces the change 
in the electron density: 



8v ext {R) = v ext exp(iq • R) + c.c. 
Sp(R) = pexp(iq ■ R) + c.c, 



Po 



(18a) 
(18b) 

(18c) 



III. STATIC LINEAR RESPONSE OF 
HOMOGENEOUS FERMIONS AND 
ADJUSTMENT TO THE LINDHARD 
STRUCTURE FACTOR 

The static properties of a homogeneous electron 
liquid are determined by the structure factor, x(o)- 



In the ground state of a homogeneous liquid the non- 
zero values at the N — 3 level are the density po and 



the second and the fourth order dynamic tensors, $ 



(2) 



Poc 2 (po)8 i3 , $ 



(4) 
ijkl 



PoC4(po)5 ■ Sijki- In a static linear 



response problem all the odd-order kinetic tensors remain 
zero. Therefore, to study the static linear response of 
the system we let the values of p, $^ 2 ^ and <E>( 4 ' vary 
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harmonically in space around their stationary values: 



where 



^ext 

P 



-_' Xt e q + C.C, 



(2) 



(4) 
ijkl 



c 2 poS lj 



pe iqR + c.c 



(■ 



+ c.c 

3c 4 po%fe + (^e^ fl 
9i + w ext )e 



+ c.c. 



iq-R 



+ C.C. 



(19a) 
(19b) 

(19c) 

(19d) 
(19e) 



where the underlined variables are the linear response 
coefficients, while 



, ^ dV xc {p) ( , 
V{q) = f~ -^p- {po) - 



(20) 



With the use of Eq. (fT6|) the infinitesimal deviation of $' 4 ^ 
has the following form: 



ZDpSySu + Qc 2 8idS + Po# •« > ( 21 



6 (4) 



D = 



<9(p(c 4 - c§)) 



9p 



_ 15*' 



Next we consider what terms can be in (j)( 4 \ Our pur- 
pose is to make sure that the static response in the non- 
interacting case would resemble the Lindhard function 
(static density response of free electrons). The terms 
added should include the derivative of the available quan- 
tities, i.e., the density and the stress tensor, so that 
they will be vanishing for uniform densities. Further, 
since is a fourth-order tensor, it needs to be con- 
structed from available tensors; the only ones available 
in the static limit are V,, p, &ki and Sji. It is easy to 
see by inspection that only the following local terms are 
available to first order in the perturbation and to lowest 
orders needed in V,-: 



■\f kl {R) = AViViVfcVip - - 6c 2 h8 l] V k Vip, (22) 

where A, / and h are dimensionless parameters. Even in linear response these terms can be augmented by terms 
involving further derivatives, e.g., terms involving a Laplacian of the components in Eq. (|22p , (i.e., q 2 in Fourier 
space) but as orbital- free methods should be primarily geared towards the long- wavelength limit, we do not consider 
here such higher order terms in q. 

The additional terms yield the following relation between the linear response coefficients p, $ < - 2 - ) , and $' 4 ' : 



Finally, the linearized equations read 
and 



iqjqkqi +6c2hpqiqj5 k i. (23) 



(2) 

Qa!ki a + QiPQVp + ^PolZcxt = °> ( 24a ) 



3(c 2 % + /<ftfc)(?a£&) + 3(c 2 + fq 2 )q t $Xl 

+ ^3 (D + c 2 p V + hc 2 q 2 ) S^qu + \jPoV + M 2 + 3hc 2 ^J qiqjqk^j p 



+ [ic 2 5 tJ q k + -qiqjqkj pav cxt = 0. (24b) 

where the index a is summed over. The only preferential direction in the problem is the momentum vector, q, so that 
the dynamic tensor, &- 2 \ can be decomposed into the two form-factors: 



$]f=^(0) + «f$W. 
3 q 2 ~ 

Upon substituting this resolution into the initial equations and equating independent spatial tensor components we 
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arrive at three equations for p, $ (0) and J> (1) 

$(°) + $(!) + po {}p = -p^ 




■cxt " 



Upon solving these linear equations one gets: 



1 



-x- 1 = 



-Xo 



where 



7T 2 1 + 20(2/- 1/12)?7 2 



fc F X ° 1 - 2Ahi] 2 - 8OA77 4 ' 
and we introduced the dimensionless momentum: 



(25a) 
(25b) 



(25c) 



V 



g 

2k F ' 



Eq. (|25bp is the definition of the structure factor renor- 
malized with respect to two-body interactions. There- 
fore Xo should be the structure factor of non-interacting 
electrons. In order to adjust our theory to the realistic 
description of electrons one should compare xo to the 
Lindhard function: 



1 1 



k F XUnd ~2 + 4t? 



hi 



1 



V 



1 



V 



The freedom in choosing the parameters A, /, and h al- 
lows us to fit our structure factor to the Lindhard func- 
tion. The Lindhard function has the following properties: 



-, — XLind 
Kp 



It 2 . 

2 : 



J] — > OO 
77 -> . 

r, = l 



(27) 



In order for our function, xo, to possess these properties 
we should choose: 



A = , f = — , and h = . 

80' 1 20' 36 



(28) 



A comparison between the resulting structure factor of 
the proposed theory with the Lindhard function and the 
structure factor provided by the 1 /9-von Weiszacker the- 
ory is given in Fig.l. 



IV. APPLICATION TO THE GROUND STATE 
PROBLEM 

We applied the $( 3 )-theory to a ground state study of 
a non-perturbative non-homogenous jellium. We chose 



a spherically symmetric infinite electron system in the 
following positive jellium background density profile: 



Po(r) 
Ap(r) 



■Ap(r), 



(29) 



i 



where we took p^ = 0.01, C — 0.9, and r = 3 and 2 
(all in a.u.). The additional non-homogeneous part of 
jellium density Ap integrates to zero so we avoid com- 
plications connected with an overall non-neutral system. 
Alternatively this system can be viewed as having con- 
stant jellium background density, p^, but with an exter- 
nal Gaussian potential: 



V b (r) 



rlCPoo ye 



which is related to Ap(r) by the Poisson equation: 
1 d d 

^ar r2 a^ (r) = 47rAp(r) - 



(26) The Dirac exchange is used here: 



V XC (R) 



(37T 3 ) 



2)1/3 



-p(R) 



1/3 



(30) 



(31) 



(32) 



No correlation energy was employed (its contribution is 
very small; it will be included in future studies). The 
simulations were performed by adiabatical turning on the 
nonhomogeneous part of the jellium positive background 
density, Ap. Initially the electron and the jellium den- 
sities are homogeneous, poo- The odd-order kinetic ten- 

(3) 

sors, Jj and ^j^, are zero and the even order tensors are 

(2) (2) 

those of a homogeneous electron liquid = Pooe-, 



and $ 



(4) _ 
ijkl 



W 
-'ijkl 



with the e-tensors given in Eq.|T 



We then propagate the set of Eqs.([7]) while the jel- 
lium density gradually changes from homogenous to the 
final po{f)', this ensures that the system remains at the 
ground-state for all times. We implemented the adiabatic 
density by setting 

Po (R, t) = (1 - g(t))p (R) + g(t) Poo , 

where g(t) is a smooth function rising from to 1; we 
chose here, quite arbitrarily, 



g(t) 



i 



1 + exp 



((V) 
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FIG. 2: The electron density profiles for a jellium model with background positive density (lower solid line) given by Eq. ()29P 
with £ = 0.9 and rg — 2, 3 (left and right graphs respectively) for Thomas- Fermi, 1/9-von Weiszacker, $ (3) HTDFT, and the 
Kohn-Sham orbital based approaches. All quantities are in a.u. 



and used 

t = 'St. 

The width parameter, r, was typically taken as 50a. u.; 
this value was more than sufficient for adiabatic conver- 
gence. po(R, t) is then used for the definition of the time- 
dependent potential, Eq.Q. 

The evolution of the system is then determined from 
the four first equations in (O, with the <&( 4 ) tensor given 
by Ea.lflT)]) and Eq.®. The 3D equations were dis- 
cretized and the derivative were evaluated by Fourier- 
transforms, as was the Coulomb integral. Grid spacings 
of 1.6 a.u. - 2 a.u. were sufficient to converge when the 
hole width parameter, r , was set at 2.0 or 3.0 a.u., re- 
spectively. A simple fixed step Runge-Kutta algorithm 
with dt=0.2 a.u. was used to evolve the equations in 
time. 

We compared the results to Thomas-Fermi, von- 
Weiszacker, and plane-wave Kohn-Sham simulations. 
The latter were done by a standard plane wave code; in- 
terestingly, we found that the grid spacing needed to con- 
verge the Kohn-Sham plane wave simulations had to be 
smaller by about 20% than those needed in the HFDFT 
code, so that they were about 1.3 and 1.6 a.u. for r =2.0 
and 3.0, respectively. The grids contained typically (20) 3 
points. 

Fig. 2 shows that HTDFT gives essentially the same 
density as the Kohn-Sham approach, while the von- 
Weiszacker and Thomas-Fermi results deviate signifi- 
cantly. Since the two-body interaction is treated the 
same in all four simulations, this proves that the hydro- 



dynamic approach yields, even for this system which is 
shifted strongly away from uniformity, the same densities 
as the essentially exact description of the kinetic energy 
in the Kohn-Sham approach. 
V. TIME-DEPENDENT LINEAR RESPONSE 
AND THE COLLECTIVE MODES 

In our previous paper (39l | we studied the ground-state 
of a homogenous electron gas at the N = 2 level, with 
the assumption that cf>( N+1 ' > is zero. Here we extend the 
studies to N = 3, with cj)^\ as given by Eas.([2"2" J) ,([2"g ]) . 
We derive the governing formulae in general, and arrive 
at analytical limits in the long wavelength limit (where 
</> 4 is not contributing), showing new kinds of excitations. 

In the ground state of a homogeneous liquid the non- 
zero values at the N = 3 level are the density po and 

(2) 

the second and the fourth order dynamic tensors, = 

PoC2(po)8ij, (frvjlt = poc i (p )6t j 6ki- To study the linear 
response of the system we let all the values in the problem 
vary harmonically around their stationary values: 

P = Po+ (pe^-q-B) + c c N t (33a) 

Ji = l ie - i[ujt - q - R) + c.c, (33b) 

- c 2Po 6 l3 + (^fe-^'-'-^+cc.) , (33c) 

$gi = igie-^-^+cc, (33d) 

W l V = qiv(q) (ipe-^-f 1 *) + c.c.) , (33e) 
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After linearizing Eqs.([7]) one gets: 

(34a) 

iv(q)p, (34b) 

(34c) 



ujp 






= Q^tJ + 








= 3(c 2 <5y + 



3 (D + c 2 po« + frc 2 g 2 ) + f ~Po« + M 2 + 3/ic 2 J qiq^qk^j £, (34d) 

where the linearized variation of <J>^ 4 ^ is taken from Eq. (|23p . and a is again summed over. This is a system of linear 
homogeneous equations, and to find its solutions we have to diagonalize it. 

_____ /Q\ 

All the variables in Eqs. (|3"4")) could be expressed in terms of p and ___y . Therefore, we can consider the equations 
on and p only without losing any solutions. In matrix form these equations read: 



q — 



Tr (* (2) Q) + Pov(q)p, (35a) 

j 2 

q 2 



^* (2) = ( C 2 + /g 2 )(* (2) +2{*( 2 ',Q})+( C2 / + / ? 2 Q)Tr($( 2 )Q) 

+ (IZ(q) + QZ'(q))p, (35b) 



where 



with dispersion 



Z(q) = D + c 2 p V + hc 2 q 2 , 
Z'{q) = 2Z+]p d(q)q 2 +Aq 4 



lu 2 = \jhk\q 2 



3hc2q 2 



and Tr means matrix trace, Zy = <5y is the 3x3 unity 
matrix, Qij = qiqj/q 2 , curly brackets denote an anti- 
commutator, and capital bold face letters refer here to 
matrices. Without loss of generality we can always as- 
sume that the wave- vector q is directed along the x-axis 
(q = (q,0,0) T ), so that 



Q 




(36) 



There are several solutions for these equations. The first 
three solutions are decoupled from the density fluctua- 
tions so that p = for all of them. They are: 



'0 1 0\ (0 V 
$ (2) =.1 0,0 
,0 0/ 10 0; 



(37a) 



which corresponds to the dispersion relation 

lu 2 = 3/5k 2 F q 2 , (37b) 

and 



'0 N 
$ (2) = | 1 
,0 1 0, 



(38a) 



(38b) 



The first two solutions, Eq. (|37ap , correspond to trans- 
verse sound as the current is given as: 



w (_)___ = q^f = |«| , | 




(39) 



Note, however, that the velocity of this transverse sound 
mode is different from the one found for the same mode 
within the N = 2 theory [39]. The third solution in 
Eq. (|38p is a new sound mode. This mode involves nei- 
ther density nor current fluctuations and corresponds to 
transverse quadrupole fluctuations of the Fermi sea. 

The next three solutions are found by representing the 
tensor <I>( 2 ) in terms of the remaining diagonal tensors (I 
and Q): 



*( 2 ) = a Q + j3(I-Q), 



(40) 



which leads, upon insertion into Eqs. ([35|) . to the following 
equations for a, (3, p: 



u> 



a - 6(c 2 + fq 2 )a+(Z(q) + Z'(q)))p. (41a) 



-tP = c 2 a + (c_ + fq 2 ) (3 + Z(q)p, 



p = a + p v(q)p. 



(41b) 
(41c) 
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FIG. 3: The elementary excitation spectrum provided by 
Quantum Hydrodynamics (QH), HTDFT $ (2) and HTDFT 
$' 3 ' theories. QH gives only a plasmon mode. HTDFT also 
gives transverse sound modes which mimic the RPA elemen- 
tary excitations in Fermi liquid. $' 3 ' HTDFT gives addi- 
tional sound modes with resi ject to $ (2) HTDFT confirming 
the conjecture made in Ref. [391] that when increasing the or- 
der of HTDFT new sound modes should appear, and they will 
gradually cover the entire continuous RPA density of states 
of Fermi liquid. 



This set of equations has complicated solutions, which, 
however, could be simplified in low-wavelength limit. In 
this limit, we can leave only the leading terms in q; in the 
effective potential it is the divergent Fourier components 
of the Coulomb potential. In the long-wavelength limit 
the system of equations has the following form: 



6 3 
1 1 1 

a(q) a(q) 







fa 


{ 




P 









(42) 



where a(q) = 4tt p / \q 2 c 2 ) , uj' 2 = iv 2 /(q 2 c 2 ), p' = 
Airpop/q 2 , and a(q) ^> 1. Dropping the terms of or- 
der a(g)~ 1 and smaller, the three eigenvalues and cor- 
responding eigenvectors are: 



U3 



k 2 F q 2 , (a,/3,p') = (1,2/3,-1); (43a) 
k 2 F q 2 , (a,/3,p') = (1,0,-1); (43b) 



u 2 P + -k 2 F q 2 , (a,/3,p') = (0,0,1), (43c) 
5 - 



where tup — Airpg is the plasmon frequency. Note that 
the first two of the three modes (|43a[ I43b[) have the 
same eigenvalues as the transverse modes in Eqs.f37bJ 
and ([38E]) . 

The total spectrum given by N = 3 HTDFT for ele- 
mentary excitations in the homogeneous electron gas is 
given in Fig. The spectrum found differs from that 
of the N = 2 approach by an additional sound mode 



with velocity ^3/5kp and by shifting the previous sound 
modes from y/3/5kp to y/l/hkp. This result confirms 
the conjecture made in Ref. [39], that with increasing N 
new sound modes should appear, and that they will grad- 
ually cover the entire continuous RPA density of states 
in the Fermi liquid. 



VI. CONCLUSIONS 

In conclusions, we have shown that HTDFT can also be 
used for time-independent studies. We have supplanted 
our previous conjecture where we assumed that the terms 
in the equation of motion hierarchy should be terminated 
with the next relevant cumulant {i.e., <j)( N+1 ^) being zero; 
instead, we now derived from fitting the linear 

response to a HEG. The resulting set of equations (given 
at the N = 3 level by Eqs. (J22J , ||S2J) , J?} , and JTB} ) is 
closed and can be propagated forward in time. 

The linear response in the static limit is fitted to the 
Lindhard function HEG for both short, intermediate and 
long wavelengths (for comparison, the 1/9 in the von- 
Weiszacker approach is obtained to fit long wavelengths, 
while a fit to long wavelengths would have required re- 
placing the 1/9 by 1 in the von-Weiszacker theory). We 
have then applied HTDFT away from equilibrium, for a 
case of a jellium density with a deep hole in the middle, 
and have shown excellent agreement with the Kohn-Sham 
results, in a case where more approximate theories such 
as Thomas Fermi and von-Weiszacker fail; this is directly 
due to the fact that their structure factor do not follow 
the Lindhard function except at low wavelengths. 

The last part of the paper dealt with time-dependent 
linear response studies at the present level, N = 3. The 
analytical studies have confirmed our previous assertion, 
that as the level of the theory increases, more and more 
transverse excitations are found. A new excitation at the 
N=3 level couples neither to the current nor to the den- 
sity. All excitations, including the new ones, lie within 
the RPA density of stats of elementary excitations in a 
Fermi liquid. 

Future work will study the applicability of the ap- 
proach to covalent chemical systems, where the direc- 
tionality of the tensors should enable a correct descrip- 
tion even at a low N, possibly as low as N = 3. Further, 
dynamic susceptibilities will be studied so that further 
terms, depending on J, etc., will be included in the 



terminating cumulant (cfr here, 



(N+l) 



in general) so that 



the theory will be valid over a wide range of frequencies 
and wavevectors. In addition, effects of magnetic fields 
are straightforward to incorporate. The basic formalism 
developed here and in forthcoming work will then enable 
the application to dynamical problems which straddle the 
transition between molecular and nanostructures. 

We note that other applications to fcrmionic systems 
can also be envisioned. For example, by replacing the 
zeroth-order HEG density matrix with a temperature- 
dependent density matrix, and fitting the coefficients of 
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the derivative terms in the cumulant to a temperature de- 
pendent Lindblad expression, we will get a temperature- 
dependent HTDFT theory which can be applied to plas- 
mas and to studies of narrow conduction bands. Simi- 
larly, applications to nuclear systems can also be envi- 
sioned. 

Other future improvements will include better meth- 
ods to solve the time-dependent HTDFT equations. One 
approach will be to include external electric fields that 
will have dipole and quadruple (or higher) components 
that will be time-dependent. The electric fields will be 
chosen, at each time-instant, to remove energy from the 



system (i.e., to reduce the trace of $^ 2 ^ plus the total 
potential). 
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